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ABSTRACT 



^ I We simulate the characteristics of the Galactic population of radio and 7-ray pulsars 

using Monte Carlo techniques. At birth, neutron stars are spatially distributed in the 
Galactic disk, with supernova-kick velocities, and randomly dispersed in age back to 
10^ years. They are evolved in the Galactic gravitational potential to the present time. 
From a radio luminosity model, the radio flux is filtered through a selected set of radio- 
survey parameters. 7-ray luminosities are assigned using the features of recent polar 
cap acceleration models invoking space-charge-limited flow, and a pulsar death valley 
further attenuates the population of radio-loud pulsars. Assuming a simple emission 
geometry with aligned radio and 7-ray beams of 1 steradian solid angle, our model 
predicts that EGRET should have seen 7 radio-loud and 1 radio-quiet, 7-ray pulsars. 
With much improved sensitivity, GLAST, on the other hand, is expected to observe 76 
radio- loud and 74 radio-quiet, 7-ray pulsars of which 7 would be identified as pulsed 
sources. We also explore the effect of magnetic field decay on the characteristics of the 
radio and 7-ray pulsar populations. Including magnetic field decay on a timescale of 5 
Myr improves agreement with the radio pulsar population and increases the predicted 



^Michigan State University, National Superconducting Cyclotron Laboratory, East Lansing, MI 48824-1116 



- 2 - 



number of GLAST detected pulsars to 90 radio-loud and 101 radio-quiet (9 pulsed) 
7-ray pulsars. The lower flux threshold allows GLAST to detect 7-ray pulsars at larger 
distances than those observed by the radio surveys used in this study. 

Subject headings: radiation mechanisms: non-thermal — magnetic fields — stars: neu- 
tron — pulsars: general — gamma rays: theory 

1. Introduction 

With the advent of the Compton Gamma Ray Observatory (CGRO), the number of 7-ray 
pulsars has grown to eight, with several additional candidates. We anticipate that many more 
pulsed sources will be added to the list with the future telescope, Gamma-Ray Large Area Space 
Telescope (GLAST) scheduled for launch in 2006. The new body of data with a much larger 
set of statistics will be essential in further constraining pulsar models. Many questions, such as 
the mechanism for radio emission and its relationship to 7-ray emission, and the location and 
geometry of the 7-ray emission currently remain unanswered. Among the known 7-ray pulsars, 
only Geminga is radio-quiet or at least radio weak (Kuzmin &; Losovsky 1997 and Malofeev & 
Malov 1997). Geminga does not emit conventional radio emission that is detectable through the 
current radio surveys. As a result, wc refer to this object as a radio-quiet, 7-ray pulsar. Of the 271 
sources listed in the Third EGRET Catalog (Hartman et al. 1999), about 170 of these 7-ray point 
sources have not been identified with sources at other wavelengths. Recently Grenier &; Perrot 
(1999) suggested that some of these unidentified sources are correlated with the Gould Belt of 
massive stars from a nearby Galactic structure consisting of an expanding disk of gas with young 
stars (< 30 million years) inclined about 20° to the Galactic plane. Gehrels et al. (2000) found 
that the flux distribution of 120 steady sources suggests two distinct groups: one having higher 
flux, hard spectra and distributed along the Galactic plane, and the other having lower flux, softer 
spectra and correlated with the Gould Belt. The soft spectra and luminosity of this second group 
of sources are significantly different than those of the known 7-ray pulsars observed by EGRET. 
Harding & Zhang (2001) suggest that some of the sources associated with the Gould Belt are indeed 
radio-quiet, off-beam 7-ray pulsars seen at large angles to the magnetic pole. 

There arc two main types of models proposed to explain pulsar high-energy emission. Polar cap 
models (Daughcrty &: Harding 1996, Sturner et al. 1995) assume that particles accelerated above 
the neutron star polar caps produce 7-rays via curvature radiation or inverse Compton scattering 
induced pair cascades in a strong magnetic field. Outer gap models (Cheng, Ho & Ruderman 1986, 
Romani 1996, Hirotani Sz Shibata 1999) assume that acceleration occurs along null charge surfaces 
in the outer magnetosphere and that 7-rays result from photon-photon pair cascades. Polar cap 
and outer gap models predict different ratios of radio-loud to radio-quiet, 7-ray pulsars, primarily 
due to the different geometry of the high-energy emission region and its location relative to the 
radio emission region, thought to originate within tens of stellar radii of the neutron star surface. 
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Outer gap models generally predict small overlap of the radio and high-energy pulsar popu- 
lations, because the high-energy and radio pulses that are visible to the same observer originate 
from opposite magnetic poles, and large numbers of radio-quiet, 7-ray pulsars, because the outer 
gap beam is much larger than the radio beam. The study of Yadigaroglu & Romani (1995), using 
the outer gap model of Romani k, Yadigaroglu (1995), found that EGRET should have detected 
about three times as many radio-quiet pulsars as radio-loud pulsars in 7 rays. Most of the pulsar 
population should be seen only in 7 rays, as Geminga-likc pulsars. The Monte-Carlo simulation of 
Cheng &: Zhang (1998) found that only about 16% of radio pulsars should be 7-ray pulsars. They 
expect 55 Geminga-like pulsars to have been detected by EGRET, about all of the unidentified 
sources in the Galactic plane, as compared to 11 radio- loud 7-ray pulsars expected. The charac- 
teristics of the 7-ray pulsar population, as seen by EGRET, in the inverse-Compton initiated polar 
cap cascade model have been studied by Sturner k, Dermer (1996a). They predict that about 75% 
of radio-selected pulsars should be 7-ray pulsars but that only 25% of 7-ray pulsars are radio-quiet. 
Although the uncertainties in these studies are large, it is apparent that the differences between 
expected popTilations of polar cap and outer gap models are larger than the uncertainties due to 
modcl-dcpcndcnt effects. 

Since there has been no statistical study of the expected high-energy pulsar populations in the 
curvature radiation-induced polar cap cascade model (e.g. Daugherty &: Harding 1996), we have 
developed a Monte-Carlo code similar to that of Sturner &; Dermer (1996a) to simulate the radio 
and 7-ray pulsar populations in the Galaxy. While we understand that the emission geometry is 
crTicial to the modeling of the 7-ray pulsar population, we present in this study a simple geometric 
model that provides insight into the number of radio-loud and radio-quiet, 7-ray pulsars detected 
by EGRET and those expected to be detected by GLAST. The simplest geometric assumption 
suggested by the polar cap model is one in which the radio and the 7-ray beams are aligned and 
equal in solid angle. With this assumption, we model the population of radio and 7-ray pulsars in the 
Galaxy. We do not attempt to fit the observations. Rather we use standard distribution functions, 
evolution techniques, radio luminosity models and recent 7-ray luminosity models suggested by an 
acceleration-cascade model of the polar cap (Zhang &: Harding 2000). The simple model presented 
in this paper predicts a ratio of radio-quiet to radio-loud, 7-ray pulsars (1/7) that is comparable 
to the one observed by EGRET (1/8). With a greater sensitivity, GLAST is expected to observe a 
ratio of 74/76. The large increase in the ratio of radio-quiet to radio-loud, 7-ray pulsars is due to 
the greater sensitivity of GLAST to detect 7-ray pulsars at greater distances from the Earth than 
the radio surveys used in this study. We plan to have forthcoming the logical extension of this work 
that will invoke more realistic emission geometries for both radio and 7-ray beams. 

2. Monte Carlo Simulation Model 

We develop a model to simulate the production of neutron stars within the Galaxy, evolving 
their trajectories, periods and period derivatives from their birth forward in time to the present. 
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We supply the radio and 7-ray characteristics to each neutron star and filter its properties through 

radio surveys and 7-ray thresholds (in and out of plane) associated with EGRET and expected for 
GLAST. These 7-ray thresholds correspond to the number of photons required for the instrument 
to identify the object as a point source. Higher thresholds would be required to obtain sufficient 
photons to identify the object independently as a pulsed source. 

2.1. Initial Pulsar Period and Magnetic Field 

In most of our simulations, the magnetic field of the pulsar is assumed to be constant through- 
out its lifetime, although we also will explore the effect of field decay. A constant field requires that 

(PP)^/^ is equal to a constant and, therefore, implies a nonzero P. Following Bhattacharya ct al. 
(1992) and Gunn & Ostriker (1970), we have assumed that the magnetic field distribution can be 
represented by a Gaussian. However, we found it necessary to include two additional Gaussians 
below the main distribution to account more fully for the inferred distribution from the period and 
period derivative of observed pulsars. The pulsar's surface magnetic field distribution is represented 
by the expression 

PB = X;Ae-(i°g^-'°g^*)'/-?, (1) 
1=1 

where the parameters take on the following values: 



Table 1 


i 


Ai 


log(i?i) 




1 


60. 


12.65 


0.45 


2 


2. 


11.9 


0.6 


3 


.001 


10.4 


4 



An array normalized to unity is created with this distribution from log{B) = 9.5 to 13.5 in steps 
of 0.02, and used to randomly select the pulsar's initial magnetic field. 

The birth rate of neutron stars is assumed to be constant during the history of the Galaxy. 
Therefore, the age of the pulsar is randomly selected from the present to 10^ years in the past. 
While most of the very old pulsars will not be observed due to their very long periods, there will be 
some observed with ages close to 10^ years as a result of their small magnetic fields and, therefore, 

small period derivatives. We use the expression by Shapiro Sz Teukolsky (1983) and Usov & Melrose 
(1995) for a uniformly magnetized neutron star where the inferred magnetic field is given by the 
star-centered dipole relation, with magnetic moment /x = BqE? /2, 

B12 = 6.4 X )V2^ (2) 

where P12 is in units of 10^^ Gauss. This relation is not used in the pulsar catalog by Taylor, 
Manchester & Lyne (1993) and in Lyne &; Graham-Smith (1998) who instead use the approximate 
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magnetic dipole moment, /j, = BqR^. Therefore, we multiplied the magnetic fields in the pulsar 
catalog by a factor of 2 in order to make comparisons. Integrating this expression over time results 
in a pulsar period in seconds at the present time by the expression: 

Pit) = [P^ + 1.54 X 10-' Bf^t] , (3) 

where t is the age of the pulsar in years and the initial period Po is in seconds. The initial period of 
the pulsar at birth is assumed to be a fixed 30 ms. The simulated distributions are rather insensitive 
to the initial period as long as Pq ^ 100 ms. Having the P{t) along with the pulsar's magnetic field, 
the pulsar's period derivative can be obtained assuming a braking index of 3 from equation (2). 
Various studies (for example, Lyne, Pritchard & Smith 1988, Kaspi et al. 1994, and Boyd et al. 
1995) measuring the second period derivative suggest a lower braking index. However, systematics 
are difficult to obtain as second derivatives of pulsar periods are very time consuming to measure. 
At this time, we have not attempted to model anything other than a braking index of 3. 

Using the recent study of Zhang, Harding &: Muslimov (2000), we have introduced a pulsar 
death line predicted by a multipole magnetic field configuration near the stellar surface within a 
space-charge-limited fiow model and described by the expression 

logP = 21ogP- 16.52. (4) 

The model includes the effect of general relativistic frame-dragging discussed in Muslimov Sz Tsygan 
(1992) and Muslimov &: Harding (1997) essential in the development of the electric field parallel to 
the magnetic field. A more recent study of Harding &: Muslimov (2001) has shown that the death 
line described by equation (4) roughly defines the boundary of pair production for a dipole field 
when non-resonant and resonant inverse Compton radiation processes are considered. Only pulsars 
with period derivatives greater than those of equation (4) are further considered. 

In addition to the death line, we have introduced a death valley in the P — P space. Otherwise, 
we find too many pulsars are predicted near the death line that are clearly not present in the 
observed distribution. In order to define the death valley, we have used a second line from Zhang, 
Harding &; Muslimov (2000) that describes the death line in the space-charged-limited flow model 
for a purely dipole magnetic field near the stellar surface given by 

log p = j log P - 14.56. (5) 

The above expression roughly describes the boundary of pair production for curvature radiation 
photons in a pure dipole field (Harding & Muslimov 2001). We have chosen an exponential decline 
in the pulsar population along constant magnetic field in the P — P plane. Wc define a distance in 
PP from the multipole death line to the location of the simulated pulsar in the death valley. We 
randomly select whether the particular pulsar is counted or rejected from further consideration in 
the code. The role of the death valley will be further discussed in a later section. 
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In addition to this case involving no magnetic field decay and a pulsar death valley, we have 
simulated a case in which the field is allowed to decay exponentially with time constant, td, given 
by the expression 

5(t) = S„i2e-*/^^, (6) 

where B012 is the magnetic field of the pulsar at birth. In the field decay case, we found that the 
pulsar death valley is no longer required. Assuming magnetic dipole spin-down and initial period, 
Po, the period and the period derivative at the present time can be obtained by 

p2 = p2 ^ y ^ 10-^52^2^13(1 - e-*/^°), and 

P = 2.44 X 10-^*5 e^*/^^-^, (7) 

where P and Pq are in seconds and t and td are in years. Due to the field decay, the spin-down 
age of the observed pulsars has to be determined using the equation 



Age = In 



3.17 X 10-^P 



Ptd 



+ 1 



(8) 



In the limit as td goes to infinity, the age becomes the traditional characteristic age of a pulsar 
given by P/2P. Due to the field decay, we find that a single Gaussian is sufficient to describe the 
initial magnetic field distribution at birth of the majority of the pulsar population in the P — P 
space. The values are listed in Table 2 for the Gaussian parameters in equation (1). 

Table 2 



i Aj \og{Bi) Gj 
1 1.0 12.75 0.4 



The role of field decay will be discussed in a later section along with the results. 



2.2. Spatial Distribution of Pulsars 

In a cylindrical coordinate system with the origin at the Galactic center, we assume that 
the birth location of neutron stars is well described by the following distributions as indicated in 
Paczyhski (1990): 

Pz{z)dz = e'^^^^^^^'p/zexpciz, and 
PR{R)dR = aRC-^/^^^R/R^^pdR, (9) 

where z is the distance from the Galactic disk and R is the distance from an axis through the 
Galactic center perpendicular to the Galactic disk and is given by: 



R^ = x'^ + y^. 



(10) 
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The constants are defined as 

an = /^-p(l + i?„ax/i?exp)]"', 

i?cxp = 4.5 kpc, 

Zcxp = 0.075 kpc, and (11) 

-Rmax = 20 kpc. 

Defined in this fashion the following integrals are normalized to unity, 

oo 

J e"l^l/'^«''p/zexpd-2; = 1, and 



-Rmax 

J aRe-^/^^-^R/Re.pdR=l. (12) 



With these distributions, pz, and pR, the initial position of the neutron stars can be chosen using 
a random number, TZ 

2; = -ln(l-7^)^exp, (13) 

with the sign of z chosen by a second random number. However, as the inversion of the integral of 
the pjt function involves solving a transcendental equation, the distance R is chosen randomly by 
creating a normalized array, 

Ij^{R) = e-^l^^^-^R, (14) 

and performing a linear interpolation between neighboring values. The azimuthal angle, ^, is chosen 
randomly between and 27r. 



2.3. Supernova Kick Velocity Distributions 

The initial velocity distribution given to a neutron star during a supernova has been the subject 
of much discussion. The three-dimensional space velocities of neutron stars described by Lyne & 
Lorimer (1994) have been obtained from a two-dimensional distribution, ^^^'^3,^ (see also Mollerach 
& Roulet (1997)) where x is proportional to the transverse velocity. The three-dimensional space 
velocity distribution can be accurately represented by an expression from Sturner k. Dermer (1996b) 
and has the form 

= ; (it?) • (^^' 

where C, = v/350 km/s and v is the random three-dimensional velocity of the neutron star in its rest 
frame distributed isotropically. The most probable and mean ( from this distribution are 0.76 and 
1.41, respectively. The advantage of this function is that it is not only normalized, but the integral 
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can be easily inverted to obtain the velocity directly from a random number with the expression 




As discussed later in the text, we have chosen this functional form, but we have used instead 
C = 'i;/120 km/s in order to obtain better agreement with the z distribution of the pulsars. 



2.4. Galactic Gravitational Potentials 



We adopt the negative of the potential functions as defined in Paczyhski (1990), 

-GMi 



and 



-GM, 



l + iln(l + f)4ln(l + g)-^tan-i 



(17) 



where $i corresponds to the spheroid, $2 corresponds to the disk potentials and is the halo 
potential. We use the same choice of parameters as in Paczynski (1990) 

ai = 0, 61 = 0.277 kpc. Mi = 1.12 x 10^° M©, 
02 = 3.7 kpc, 62 = 0.20 kpc, M2 = 8.07 x lO^^Mo, (18) 
rc = 6.0 kpc, Mc = 5.0 x IO^^Mq. 

where the variable, r, corresponds to the radial distance from the Galactic center: 

r = Vx2 + y2 + ^2_ (19) 

We have added a constant term to the halo potential that has been left out of the expression in 
Paczyhski (1990). The constant negative term dominates over the r dependent term resulting in an 
overall negative potential. In the constant term, R^, is the radius of the halo with a typical value 
of 41 kpc (Binney & Tremaine 1987). However, this term does not affect the pulsar trajectories 
as the equations of motion depend on the derivatives of the potential. It is important, if one is 
interested in the total energy to determine if the stars are unbound. Of course, if the star escapes 
the halo, the potential would change. The pulsars of interest in this study are well within the halo, 
and, therefore, this form of the potential is appropriate to evolve the trajectories. The Lagrangian 
in units of energy per unit mass has the form 

1 



L = T-U 



^) -$(i?,z). 



(20) 



where $(i?, z) is the total potential energy per unit mass. Since the Lagrangian is independent of 
0, the angular momentum, along the z direction is a constant of motion. 



V, 



cir 



so that 



(21) 
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The circular velocity, Vdr, is the motion of the star in the Galactic plane (for z 
of motion to be integrated for the R and z directions have the forms 

z = -— 

dz ' 

Of coTirse the supernova explosion will impart an initial random velocity (eqn [16]) that needs 
to be taken into account in determining the initial angular momentum and the total energy of 
the system in its Galactic orbit. Once the initial conditions are established, the trajectories are 
integrated using a fourth order Runge-Kutta routine designed to maintain a high level of accuracy 
in the conservation of total energy of one part per 10^ during the trajectory from its birth to the 
present time. 



= 0) . The equations 



(22) 



2.5. Radio Luminosity 

We obtain the radio luminosity of the pulsars at 400 MHz, L400, using the radio pulsar model of 
Narayan Sz Ostriker (1990) where the normalized luminosity distribution is dithered by a function 
given by 

PL,oo = 0.5X^e-\ (23) 

where 

A = 3.6 f log ( + 1.8] , (24) 



V (-^400) / 

and where the average radio luminosity, (-L400), is given by 

log(L4oo) =6.64+^ log (25) 

This is the luminosity law that was adopted by Bhattacharya et al. (1992) along with the above 
dithering function. This luminosity function derived from P and P was obtained by Proszyhski &i 
Przybicieh (1984). We build an array with the distribution indicated by equation (23) from A = 
to 20 in steps of 0.02, and using a random number, we linearly interpolate between neighboring 
array elements to select a A. Given the period and period derivative, equation (25) gives the 
average luminosity at 400 MHz, (L400)! and together with equation (24) provides the luminosity 
of the pulsar at 400 MHz, L4oo(mJy • kpc^). However, as discussed later, we find that we had to 
modify the intercept of equation (25) in order to achieve better agreement with the observed pulsar 
distributions. The flux in mJy is then obtained from 

where d is the distance from the Earth in kpc. We have assumed a constant solid angle of ACl = 1 
steradian. The simulated flux at 400 MHz is scaled to the observing frequency of the surveys we 
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model using a spectral index of —1.7. This is the average spectral index in the frequency range 
between the fluxes of 5*400 Siud S'1400 of the select group of observed pulsars and is in agreement with 
Johnston et al. (1992). Pulsars with fluxes greater than the survey flux threshold are detected. 



2.6. Radio Detection 

The sensitivity of a particular survey is a function of several parameters usually given by the 
expression (Dewey et al. 1985) for S'min ipiJy)., 

cr _ Cthres [Tree + Tg^yjl, b)] j W . . 

where Cthres is the detection threshold S/N, T,.ee(^) is the receiver noise temperature, Tsky{l, b)(K) 
is the sky temperature in the direction being searched, G (K/Jy) is the telescope gain, which we 
adjust to reflect various system losses, Np is the number of polarizations (usually 2), B (MHz) is 
the total bandwidth, t (s) is the integration time, P (ms) is the pulsar period and W (ms) is the 
effective pulse width. 

Expressions for S'min in the literature sometimes (Johnson et al. 1992, Lyne et al. 1998, 
Manchester et al. 1996) include an additional term, /?, in the numerator to explicitly account for 

the losses due to digitization and other system losses. In such cases, we incorporate the reported 
values of (3 by adjusting the value of G (Table 4), increasing S'min by effectively reducing the gain. 
We use the sky temperature at 408 MHz determined by Haslam et al. (1982) that has been made 
available in machine readable form (http://skyview.gsfc.nasa.gov). Given the Galactic longitude 
and latitude, we interpolate a 512 x 1024 table of sky temperatures at 408 MHz. We then scale 
the temperature to the observing frequency u (MHz) through a power law given by (Johnson et al. 
1992) in the form 

408 MHz^ 



Tskyi^) = Tsky, 408 y ' j ■ (28) 

The effective pulse width W is given by 

2 2 2 2 2 / \ 

W =Wg + Tg^j^p + TDM + Tscat + TtrailDM, (29) 

where Wo is the intrinsic pulse width (here assumed to be 0.05P(ms)), Tgamp represents low pass 
filter time-constant applied before sampling or, if unknown, an assumed value at twice the sampling 
interval (2 At), tum is the dispersion smearing t imc over one frequency interval, Au (MHz), Tg^at 
(ms) is the time broadening of the pulse due to interstellar scattering. The dispersion broadening 
time, tbm (ms), across one frequency channel, Az/ is related to the dispersion measure (DM) and 
has the form ^ ^ , ^ 

TDM = 8.3 X 10^ 5 — . (30) 

The dispersion measure DM(pc/cm^) is obtained using the Taylor h Cordes (1993) distance model, 
where we have translated the FORTRAN routine into C code and used an extended trapezoidal 
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integration routine (Press et al. 1992). The distance model also provides the scattering measure 
SM (kpc • m"^"^/^), which allows one to obtain the broadening time due to interstellar scattering, 
Tscat (ms), used in equation (29) with the expression 

w = 1000(— ) «(— ) , (31) 

where the scattering time is scaled (Johnston et al. 1992) from 1 GHz assumed in the Taylor & 
Cordes (1993) model to the observing frequency, v. The final term in (29), TtrailDM, corresponding 
to the fourth term in the expression for W in Dewey et al. (1985), is the additional time broadening 
when the sampling is performed at a DM different than the actual DM of the pulsar. This term 
becomes important for low period pulsars. Since we are interested in pulsars with periods greater 
than 30 ms, we have neglected this term in our simulations. As will be indicated later, we do not see 
many pulsars with periods less than 100 ms where this term might add an important contribution 
to the smearing time, an effect that we will consider including in subsequent refinements of the 
model presented here. 

We have selected a set of eight surveys to filter our simulations in making compar- 
isons with the 707 pulsars of the Princeton catalog (Taylor, Manchester & Lyne 1993) from 
http://pulsar.princeton.edu/pulsar/catalog.shtml. We selected pulsars from the catalog with posi- 
tive period derivatives, periods greater than 30 ms and pulsars that are not in globular clusters and 
are not in binary systems to form a comparison set of 445 pulsars observed by the indicated surveys 
in Table 3. This group includes 90% of the 496 pulsars having these characteristics and observed 
by all the surveys in the Princeton catalog. We have multiplied by a factor of 2 the magnetic field 
strengths of the pulsars in the catalog to compare with our adopted surface field in equation (2). 

Table 3 



Selected Pulsar Surveys 



Survey 


Name 


Reference 


Eff. 




Boundaries 












Type 


a/£ 


S/b 


1 


Molonglo 2 


Manchester et al. (1978) 


1.00 


Eq. 


[0,360] [ 


-85,20] 


2 


Green Bank 2 


Dewey et al. (1985) 


0.22 


Eq. 


[0,360] [ 


-18,90] 


3 


Green Bank 3 


Stokes et al. (1986) 


0.17 


Eq. 


[0,360] [ 


-18,90] 


4 


Arecibo 2 


Stokes ct al. (1986) 


0.58 


Gal. 


[40,65] [ 


-10,10] 


5 


Parkes 1 


Johnson et al. (1992) 


0.91 


Gal. 


[270,20] 


[-4,4] 


6 


Arecibo 3 


Nice et al. (1993) 


0.49 


Gal. 


[35,65] [- 


8.0,8.0] 


7 


Parkes 2 


Manchester et al. (1996) 


0.95 


Eq. 


[0,360] 


[-90,0] 






Lyne et al. (1998) 










8 


Jodrell Bank 2 


Clifton k Lyne (1986) 


0.93 


Gal. 


[-5,105] 


[-1,1] 






Clifton et al. (1992) 











The essential parameters (Z. Arzoumanian, priv. comm.) of these surveys are indicated in Table 
4 and discussed in the text. The gain for the Jodrell Bank 2 (A. Lyne, priv. comm.) has been 
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adjusted to reproduce the long-period 5min in Clifton et al. (1992). We also do not adjust the gain 
as a function of zenith angle for the Arecibo surveys. 

Table 4 



Survey Parameters 



Survey 


Gain 


Crhres 


T 

rec 


u 


t 


At 


r 


B 


Ai/ 




(K/Jy) 




(K) 


(MHz) 


(s) 


(ms) 


(ms) 


(MHz) 


(MHz) 


1 


5.10 


5.4 


210 


408 


40.96 


20 


40 


3.2 


0.8 


2 


0.89 


7.5 


30 


390 


137 


16.7 


33.5 


16. 


2.00 


3 


0.95 


8.0 


30 


390 


131 


2.0 


2.2 


8.0 


0.25 


4 


10.9 


8.0 


90 


430 


39.3 


0.30 


0.4 


0.96 


0.06 


5 


0.24 


8.0 


45 


1520 


157 


1.20 


2.4 


320 


5.0 


6 


13.35 


8.5 


75 


430 


68.2 


0.52 


0.5 


10.0 


0.078 


7 


0.43 


8.0 


50 


436 


157 


0.30 


0.6 


32.0 


0.125 


8 


0.40 


6.0 


40 


1400 


524 


2.00 


4.0 


40 


5 



These parameters are used to calculate a minimum radio flux from equation (27) that is compared 
to the calculated radio flux from equation (26) scaled to the observing frequency with a spectral 
index of —1.7. The efficiency of the survey (Eff. in Table 3) has been obtained from the number of 
survey beams times the reported solid angle at half-power beamwidth, divided by the area enclosed 
by the survey boundaries on the sky. The efficiency is assumed to be constant over the area 
surveyed. This is not the case for some of the surveys where more sampling was done in certain 
regions while other regions were sampled more sparely. In Figure 1, we indicate the flux thresholds 
for the different surveys used in the simulation assuming a sky temperature of 200 K at 408 MHz 
characteristic of the Galactic disk and a dispersion measure of 200 cm~^ • pc. We realize that for 
some surveys that do not cover this Galactic disk, this temperature might not be appropriate. We 
choose a common temperature for the sake of comparing the Smin of each of the surveys at the 
corresponding observing frequency. However, in the actual simulation the sky temperature at 408 
MHz is obtained using the all-sky map of Haslam et al. (1982) and, then scaled to the particular 
observing frequency of the survey being tested. For the calculation of S'min in Figure 1 only, we 
assume a Tgcat obtained from the formula by Bhattacharya et al. (1992) scaled to the observing 
frequency and given by 



Tscat 



2Q-4.62+1.141og(£)M) _j_ -|^Q-9.22+4.461og(DM) 



/400^ 



In the actual simulation, we obtain the scattering measure, 5M, from the Taylor & Cordes (1993) 
distance model as indicated above. In Figure 1, The smallest flux thresholds occur for large periods 
where the intrinsic pulsar width. Wo, dominates over other terms in equation (29) resulting in a 
constant S'min from equation (27). The flux threshold increases with decreasing period when the 
intrinsic pulse width is dominated by the other pulse broadening terms in equation (29) and also 
as a result of the period dependence in equation (27). The Green Bank 2 and Molonglo 2 surveys 
are not very sensitive to pulsars with periods smaller than 0.1 seconds. 
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2.7. j-ray Luminosity 

We have taken the expressions describing the 7-ray luminosity from the work of Zhang & Hard- 
ing (2000) where a polar-cap model simulates the pair cascade region, with curvature radiation of 
the primary particles and synchrotron radiation and inverse Compton scattering of subsequent 
higher generation pairs. In addition, the model uses the self-consistent acceleration model of Hard- 
ing & Muslimov (1998) to produce the primary particles. The 7-ray luminosity is constrained by 
equation (61) of Zhang & Harding (2000) to be less than the spin-down luminosity, or 

-Bp,i2-P~^^^(cosa )-5/4 > 1.65, (33) 

where we have used a = 30°. If equation (33) is not satisfied, the 7-ray luminosity is equal to the 
spin-down luminosity, 

L,a = 9.68 X 10^''Bl-,^p-%5, (34) 

where I45 is the moment of inertia in lO^^g • cm^. If equation (33) is satisfied, then the pulsar is in 
regime I if 

^1/7 p-9/28(^^g ^)3/28 > g (35) 

otherwise the pulsar is in regime II. The 7-ray luminosities for the two regions are then given by 
the forms 

L^{I) = 4.8 X 10^^ erg • s"^ B^^J^P^^'^/^^icos af^^ , and 

L^(//) = 1.6 X 10^^ erg • s"^ 5p_i2^"^^^cos a)^/^. (36) 



These expressions correspond to the equations (59) and (60) in the work by Zhang &; Harding 
(2000) with an improved form for the luminosity in region I based on a more accurate expression 
for I . Equation (36) reproduces the dependence of 7-ray luminosity on E^^'^ oc BpP"'^ observed 
for CGRO detected pulsars (Thompson ct al. 1997), and also agrees with the empirical luminosity 
law for CGRO pulsars found by McLaughlin & Cordes (2000). Having the 7-ray luminosity, we 
can compute the flux assuming a solid angle AO, by the form 

= A^' ^^^^ 

where we have assumed the same solid angle of 1 steradian as in the case for the radio emission. 
We have used the flux thresholds (D. Thompson, priv. comm.) for in and out of the 
Galactic plane as indicated in Table 5. 

Table 5 

In and Out-of-Plane Flux Thresholds 

EGRET GLAST 

In-plane < 10" 1.6 x 10"''' photons/(cm^ • s) 5.0 x 10~^ photons/(cm^ • s) 
Out-of-plane ||6|| > 10° 7.0 x 10"^ photons/(cm2 • s) 2.0 x 10"^ photons/(cm2 • s) 
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2.8. 7-''^y Spectra 

The number of simulated 7-ray pulsars depends crucially on the average 7-ray energy assumed 
in the calculation that converts the luminosity from erg/s to a flux (photons ■ cm~^ • s~^) which 
gets compared to the instrument flux thresholds. If the 7-ray spectrum is assumed to be a power 
law of the form 

where A is a normalization factor and a is the spectral index, the average energy above an energy 
threshold, eth, can be obtained by 

(^7) = ^ (W-*.) (39) 

where e^ax is the high energy cutofi' in the spectrum. For example, typical 7-ray spectra observed 
by EGRET have a threshold of 100 MeV with spectral indices varying from 1.2 to 2.2 and the 
cutoffs vary from 5 to 30 GeV. For this range, we indicate the average 7-ray energies in Figure 2a 
as a function of both the spectral index and the cutoff energy. 

The contours of average energies in Figure 2a reflect a flat distribution between 400 to 800 
MeV where presumably most of the 7-ray pulsars exist. Rather than assuming a constant average 
energy for all simulated pulsars, we obtain a rough estimate of how the average energy varies 
with the period derivative and period of the pulsar. The spectral index of the observed 7-ray 
spectrum can be estimated from the polar cap cascade model (Harding & Daugherty 1999), which 
suggests a dependence upon the parameter P/B12. A crude estimate of this interdependence can 
be represented by an expression 

a = +0.85 - 0.45 log • (40) 

The maximum energy in the 7-ray spectrum is attributed to the photon attenuation in the creation 
of electron-positron pairs as the photon propagates through the magnetosphere. The pair creation 
cutoff energies, emax, can be derived from cascade codes. Baring &; Harding (2000) have formulated 
the dependence of the maximum energies of photon escaping the magnetosphere due to pair creation 
at different altitudes of emission. A rough representation of these energies can be made by the 
following equation 

w(GeV) = (41) 

Having the period derivative and period dependence of the spectral index and cutoff energy from 
equations (40) and (41), we are able to estimate the average energies from equation (39) as a 
function of the period derivative and period shown in Figure 2b. 
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2.9. Emission Geometry 

In this study, we have used a simple geometry for the radio and 7-ray emission assuming in 
both cases a constant emission solid angle of 1 steradian in which the radio and 7-ray beams are 
aligned and pointing in the direction of the Earth. There would be 47r other pulsars under the same 
circumstances that would not be detected because the emission cone would not be in our direction 
so we do not follow these in our code. Comparisons are, therefore, made between simulated radio 
pulsars and pulsars detected by EGRET and GLAST, observed as point sources without reference 
to the photon flux required for the identification of 7-ray pulsations. 

2.10. Implementation 

An event is initiated by selecting a neutron star with a random age, magnetic field and a 
fixed initial period of 30 ms. The period is evolved to the present time assuming dipole spin-down 
of either a constant magnetic field or an exponentially decaying magnetic field. Although the 
total number of neutron stars are counted, only events to the left of the multipole death line in 
the period derivative - period diagram are further followed. The neutron star is imparted a kick 
velocity due to the supernova explosion, and its trajectory is integrated from its birth location 
forward in time to the present. If the star is within 30 kpc of the Earth, in the viewing region of 
one of the surveys, and a random number is less than the geometric efficiency of the survey, the 
simulated pulsar is assigned a radio flux that is scaled to the observing frequency of the survey 
as well as a 7-ray flux above 100 MeV. We assume that the radio and 7-ray beams are aligned 
with an effective solid angle of 1 steradian and that the Earth is within the line of sight of the 
emission beams. Considering that there are a factor of Air other pulsars whose beams are not in 
the direction of the Earth, we can estimate a neutron star birth rate in the Galaxy. If the scaled 
radio flux is above the minimum flux (Smin) of the survey, the neutron star is flagged as a radio 
pulsar "observed" by that survey. Although it may be that several surveys "observed" the same 
pulsar, we count it as one event. Independently from the radio surveys, if the 7-ray flux is above the 
flux thresholds of EGRET or GLAST, it is flagged as a 7-ray pulsar observed by either EGRET or 
GLAST. In order to obtain smooth distributions, we continue the process until 10,000 radio pulsars 
are observed. As EGRET viewed (ftp://gamma.gsfc.nasa.gov/pub/PULSAR) each of the observed 
radio pulsars in our selected group, we normalize the distributions to the number of pulsars in the 
select group of observed radio pulsars. Since there are 445 observed radio pulsars in our select 
group that represents 90% of all the pulsars in the Princeton catalog with the features that we have 
specified, this normalization allows for the determination of neutron star birth rates. We compare 
the number of simulated 7-ray pulsars predicted to be detected by EGRET with the number that 
EGRET actually detected, and also predict the number of 7-ray pulsars that GLAST should be 
able to detect as point and as pulsed sources. In this study, we define radio-quiet 7-ray pulsars 
as those which are detectable by 7-ray telescopes but not by the radio surveys. In this manner, 
we are able to tag the events as radio-loud and radio-quiet, 7-ray pulsars "observed" by EGRET 
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and/or GLAST. Of particular interest is the number of predicted radio-loud and radio quiet, 7-ray 
pulsars, how they compare to those observed by EGRET as point sources and where each of these 
source groups are located in the P — P diagram. 

3. Results 

For the purpose of this study, we have focused on performing simulations for two different cases: 
one in which the magnetic field remains constant and a pulsar death valley in the period derivative- 
period diagram is assumed and one in which the magnetic field is allowed to decay exponentially 
without a death valley. We have placed greater emphasis on the case where the field is constant. 

In Figure 3, we present a comparison of the observed (3a) and simulated (3b) distributions 
of 445 pulsars shown in Galactic coordinates as Aitoff projections. The simulation assumes no 
field decay. The distributions are governed in part by the survey regions that have been chosen, 
but also by the many primary distributions such as birth location and Galactic potentials, etc. 
As can be seen, the calculated and observed distributions are very similar. In this preliminary 
study, we do not expect to reproduce the exact numbers observed by each survey. The number 
of pulsars observed depends crucially on the assumptions made in the determination of Smin and 
the geometric observing efficiency. As mentioned, we have neglected the last time-smearing term 
in equation (29), we have assumed a constant duty fraction of 0.05 for the intrinsic pulse width 
for all periods and all surveys, we have assumed a constant receiver temperature, and we have 
assumed that the efficiency is constant over sky boundaries of the surveys. We hope to improve 
these assumptions in a subsequent study that includes a realistic representation of the emission 
geometries for both the radio and 7-ray beams. Yet we affirm that this preliminary investigation 
makes important contributions to the understanding of pulsar emission and paves the way for 
further studies. 

In Figure 4, we compare the period-derivative versus period (P — P) plot of the select group 
of observed pulsars (4a) with the one of simulated pulsars (4b). The dotted lines are shown for 
the locus of constant magnetic field with the indicated strength. The dashed lines represent the 
indicated ages of pulsars assuming a dipole spin-down field. The solid lines show the pulsar death 
lines for dipole and multipole magnetic field distributions in the curvature-initiated, space-charge- 
limited-flow model (SCLF) (Zhang, Harding & Muslimov 2000). The radio pulsars observed (4a) 
and those simulated and filtered (4b) through the select group of surveys are indicated with solid 
dots. 

Radio-loud and radio-quiet, 7-ray pulsars are represented with solid triangles and open circles, 
respectively. The observed radio-quiet, 7-ray pulsar is Geminga and a list of radio-loud, 7-ray 
pulsars detected by EGRET includes Vela, the Crab, B1951-F32, B1706-44, B1509-58, B1055-52, 
B0656+14 and J1048-5832. The pulsar J0218+4232 with a period of 2.32 ms and a period derivative 
of 7.5 X 10~^° s • s~^ is outside the bounds of the plot and of the characteristics of the select group 
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of pulsars. The simulations for this case result in 7 radio-loud and 1 radio-quiet, 7-ray pulsars 
"observed" by EGRET (shown in Figure 4b) and 76 radio-loud and 74 radio-quiet, 7-ray pulsars 
"observed" by GLAST. The previously mentioned death valley operates between these two lines 
where the code exponentially (oc e*^"'^/'^)) removes excess pulsars with a decay constant of r = 1.0 
in a random fashion depending on the pulsar's distance, d, from the multipole death line along a 
constant magnetic field. Figure 5 shows the simulated distribution using the same set of parameters 
as those that generated Figure 4b except without the death valley in effect. 

Clearly the excess of pulsars near the multipole death line in Figure 5 is not seen in the 
observed distribution in Figure 4a. However, there are still significant differences even with the 
presence of the death valley in comparing Figures 4a and 4b. The observed distribution of radio 
pulsars is narrower in period, elongated up and down along the period derivative axis, whereas the 
simulated distribution is elongated left and right along the period axis. The observed distribution 
has more pulsars with smaller period derivatives and fewer high field pulsars in the death valley 
than simulated distribution. The observed distribution of low field pulsars is narrower in period 
than that simulated, while the distribution of high field pulsars is broader than that simulated. 
Broadening the primary high field Gaussian (Table 1) or including a higher field component only 
increases the number of high field pulsars bunched up near the multipole death line and cannot 
reproduce the group of high field pulsars observed with periods between 0.1 and 2 seconds in Figure 
4a. 

The observed group of 7-ray pulsars is on average younger, with larger P and higher fields, than 
the group of simulated pulsars. However, the one observed radio-quiet, 7-ray pulsar is older with 
a higher field than the 1 predicted. An important factor in the selection of the 7-ray pulsars is the 
assumed flux thresholds (Table 5) . The comparison of the observed and simulated distributions in 
Figure 4 suggests that there are fewer observed pulsars with higher periods than in the simulated 
distribution. It might be that in addition to the minimum radio flux dependence on the period 
described in equation (27) there is a high period limit to the radio surveys where typically the 
radio sensitivity increases with increasing period. However, there is in some published figures a 
slight increase in Smin for pulsars with periods greater than a few seconds (Dewey et al. 1985 and 
Stokes et al. 1986), perhaps due to computations of the actual telescope Smin during the observing 
period. As there are no details given in these publications, we have not taken this into account in 
our model simulations. 

The clear absence of high-field, high-period observed pulsars in Figure 4a indicated by the 
shaded circle labeled HB (High field) as compared to those simulated in Figure 4b, might be 
suggestive of the decay of the magnetic field. There are observed pulsars in the region indicated 
by the arrow that should have led to pulsars in the HB region. In fact, the whole wedge-shaped 
distribution of observed pulsars in Figure 4a could be explained by field decay, including the absence 
of low-field observed pulsars in the shaded region labeled LB (low field) . If the field does not decay, 
one would expect pulsars in this region (LB) to have led to the observed low-field pulsars indicated 
by the arrow. The excess of simulated high-field, high-period pulsars have ages of the order of 
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10 years. Therefore, a decay constant of this order is required for these high-field pulsars to have 
their fields decay by an order of magnitude. In Figure 6, we present a simulation in which we 
have included the decay of the magnetic field with a time constant equal to 5 x lO*' years. In 
this simulation, we removed the death valley to facilitate comparison with Figure 5 where we have 
simulated pulsars without the pulsar death valley to directly see the effects of the field decay. With 
field decay, the pulsars do not seem to bunch-up along the multipole death line as in Figure 5 
as the field decreases with age, producing a wedge-shaped distribution at low fields, small period 
derivatives and mid-periods. 

We found that the multiple low-field and high-field Gaussians used in the main simulations 
(see Table 1) were not necessary to achieve a distribution comparable to the observed distribution 
in Figure 4a. The improvement is quite noticeable. The wedge-shaped distribution of the observed 
pulsar distribution is explained, even in the region of low period and moderate period derivative, 
where the simulation in Figure 4b has an over abundance of pulsars (LB) in comparison to the 
observed distribution. The simulations for this case result in 9 radio-loud and 2 radio-quiet, 7- 
ray pulsars "observed" by EGRET (shown in Figure 6) and 90 radio-loud and 101 radio-quiet, 
7-ray pulsars "observed" by GLAST. Introducing the pulsar death valley along with the field decay 
would remove pulsars from the death valley and yield more pulsars to the left of the valley, thereby, 
increasing the number of radio-loud, 7-ray pulsars. The simulation with field decay suggests that 
the death valley might not be a required artifact to have reasonable agreement with the number of 
radio-loud, and radio-quiet, 7-ray pulsars detected by EGRET. The main difficulty in accounting 
for the observed shape of the pulsar distribTition by field decay is being able to justify the small time 
constant of 5 x 10^ years. Several mechanisms for field decay in neutron stars have been studied 
(e.g. Goldreich Sz Reisenegger 1992). Ohmic diffusion, with a decay tinicscale ~ 2 x 10^^ yrLg/T|, 
where L5 is the characteristic length scale of the fiux loops in units of 10^ cm and Tg is the core 
temperature in units of 10^ K, dominates in fields B < 10^^ G. Hall drift in the crust, with a 
timescale ~ 5 x 10^ L\T^ / B12, dominates in fields B ~ 10^^ — 10^^ G. Ambipolar diffusion, 
with a decay timescale ~ 3 x 10^ yr L^T^ / 3^2, dominates at the highest fields. Thus, for fields 
B ~ 10^ — 10^'^ characteristic of the radio pulsar population, the decay timescale is longer than 
3 x 10^ yr. However, if the birth distribution of fields has a mean of around 2 x 10^^ G, then the 
average initial decay timescale would be 5 x 10^ yr. But since the decay timescale is a strong 
function of field strength, this initially short decay timescale would not persist as the field decreases. 
Clearly, a more detailed study is required to fully address the effect of realistic field decay in the 
radio pulsar population studies. 

In Figure 7, we compare various distributions of the indicated features of observed pulsars 
(shaded) and of simulated pulsars. The smooth simulated histograms have been obtained from a 
group of 10,000 pulsars and then normalized to a total of 445 pulsars, the number of the selected 
group of observed pulsars. The dark histograms represent the case in which there is no field decay 
and a pulsar death valley is assumed as discussed previously. The light histograms result from the 
simulation of the field decay case with no death valley. Under the assumption of field decay, the 
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pulsar age depends on the decay constant of the magnetic field (5 Myr). As a result, we show 
separate figures for the pulsar age distributions. As noted for the case of no field decay, too many 
older pulsars with long periods and with small period derivatives are produced in our simulation. 
However, the simulated distributions of radio flux and distance from Earth along with the magnetic 
field agree very well with those observed. The distributions for the case of field decay overall appear 
to better describe the observed distributions without the necessity of introducing a pulsar death 
valley between the dipole and multipole death lines. 

The period and period derivative dependence of the radio luminosity is another important 
function that crucially determines the radio selection of simulated pulsars. As mentioned previously, 
we have used the functional form developed by Bhattachary et al. (1992) who followed the work 
of Narayan & Ostriker (1990), Vivekanand & Narayan (1981) and Proszyhski & Przybycien (1984) 
and found the best fit for a luminosity function described by equation (25). While we do not 
wish to develop a fitting model in this study, we found that it was necessary to slightly adjust the 
intercept of their functional form from 6.64 to 7.2. We present in Figure 8 a comparison of the radio 
luminosity distribution at 400 MHz (log]^o(-^40o)) as a function of log^Q{P / P^) and the population 
distributions of the radio luminosity for observed and simulated pulsars displayed as histograms. 
In the logio(-^40o) versus log]^Q(P/P^) plots, we performed linear fits in order to make comparisons 
of the distributions. The resulting fit parameters are indicated in Table 6. 



Table 6 


Linear Fit Parameters 




Intercept 


Slope 


Observed 


5.2 ±0.4 


0.21 ±0.03 


6.64 intercept in Eqn. 22 


5.0 ±0.3 


0.23 ±0.02 


7.2 intercept in Eqn. 22 


5.2 ±0.3 


0.21 ±0.02 



Note that we do not expect to obtain the same parameters as in equation (25), which represents 
the average radio luminosity of the primary pulsars prior to being observed. The modified 
luminosity law with an intercept of 7.2 also fits well the luminosity distributions for 
the case assuming the decay of the magnetic field. The actual radio luminosity at 400 MHz 
is dithered about the average using the function described in equation (23) and (24). Hartman et al. 
(1997) describe how the parameters of the dithering function affect the net luminosity distribution. 
The work of Narayan &; Ostriker (1990) was done with a sample of 265 pulsars known at the 
time with measured period and period derivative obtaining a functional form that has been used 
quite often in the literature. However, it seems clear that with a much larger sample of pulsars, a 
better functional form needs to be obtained. This functional form of the radio luminosity plays an 
important role in the selection of a radio pulsar and in the ratio of radio-quiet to radio-loud, 7-ray 
pulsars. Clearly the required dithering about the average luminosity masks some dependence that 
is currently not understood, perhaps reflecting some dependence of the emission geometry on the 
period as suggested by Kijak & Gil (1998) and Rankin (1993). 
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Another important parameterization of the primary pulsars is the initial velocity kicks given 
to pulsars at the time of their birth from asymmetric supernova explosions. These velocity distri- 
butions have been discussed extensively by Dewey & Cordes (1987), Lyne & Lorimer (1994), Bailes 
(1989), van den Heuvel k. van Paradijs (1997), and Lorimer, Bailes & Harrison (1997) indicating 
space birth velocities of 400 - 500 km/s. Hansen & Phinney (1997) indicate that these studies did 
not include the selection effects from flux limits of the pulsar surveys or the accuracy of the proper 
motion determinations and conclude that a Maxwellian distribution with a mean velocity of 300 
km/s adequately describes the observations. More recently, a study undertaken by Cordes & Cher- 
noff (1998) concludes that the three-dimensional velocity distributions can be best accounted for 
using a two component Gaussian function with mean velocities of 175 and 700 km/s representing 
86 % and 14 % of the population. In this study, we have chosen to use the functional form of Lyne 
& Lorimer (1994) with some modification. We used the z distribution of the radio pulsars from the 
Galactic disk to establish the kick velocity parameter. In Figure 9, we present the z distribution 
(shaded) of observed pulsars from our select group of 445. 

The spikes at « 1.7 kpc reflect the limitation of the distance model (Taylor & Cordes 1993). 
The thin and thick solid histograms represent the simulated distributions using C, = t^/120 km/s 
used in this study and Q = u/SSO km/s used by Lyne & Lorimer (1994), respectively. Since 
the primary z distribution of pulsars is assumed to be exponential, we have fit the observed and 
simulated distributions shown in Figure 9 with exponential forms obtaining the following widths 
shown in Table 7. 



Table 7 

Widths (kpc) of Exponential Fits of the z Distribution 
Observed For C = f/120 km/s For C, = i;/350 km/s 
0.18 ±0.01 0.18 ±0.01 0.34 ±0.01 



The functional form of Lyne & Lorimer with a parameter of C = 11/350 km/s results in a 
distribution with most probable and average velocities of 266 and 494 km/s, respectively, and 
leads to a significantly broader distribution twice the width of the observed distribution. Much 
better agreement is obtained with a parameter of C = f/120 km/s with a most probable and 
average velocities of 91 and 168 km/s, respectively. We do find that the z distribution is slightly 
sensitive to other parameters of our model described in this study. Again we are trying to use the 
best available distributions of pulsar features, but comparing the resulting distributions "by eye" 
has required us to make some small changes. The agreement is equally good for the case 
assuming the decay of the magnetic field with a slight narrower distribution with a 
width of 0.14 ± 0.01 kpc. As a result, we have used this same velocity distribution for 
both the cases where the field is constant and is allowed to decay. 

In Figure 10, we present a set of P — P diagrams for radio-quiet, 7-ray pulsars as seen by 
EGRET (10a) and GLAST (10b) and for radio loud, 7-ray pulsars as seen by EGRET (10c) and 
GLAST (lOd) represented by solid circles. We have also included the group of 7-ray pulsars de- 
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tected by EGRET in Figures 10a and 10c shown by solid triangles. In order to obtain smoother 
distributions of these different groups in the PP diagram, we have simulated a group of 10,000 
pulsars and depict their distributions by the contoured regions. The gray scaled background repre- 
sents the energy- integrated 7-ray and average radio luminosities in the upper panels (10a and 10b) 
and in the lower panels (10c and lOd), respectively. The average radio luminosity is proportional 
to pV3p-i^ while the 7-ray luminosity in regime II drops faster and is proportional to P^/^P 
The most intense luminosities for both radio and 7 rays arc in the upper left in the P — P diagram, 
but here the pulsars are very young and few in number. The numbers indicated in the legend 
represent the number of simulated pulsars from the group of 445 radio pulsars. The model predicts 
that GLAST should observe 76 radio-loud, 7-ray pulsars compared to 7 predicted for EGRET 
detected as point sources. The model predicts that GLAST should observe 74 radio-quiet, 7-ray 
pulsars compared to 1 predicted for EGRET. The GLAST sensitivity for blind period searches is 
expected to be about the same as the EGRET point-source detection sensitivity (S. Ritz, priv. 
comm.). GLAST will, therefore, be expected to identify 7 of these 74 objects as pulsed sources. 
The distribution of radio-loud, 7-ray pulsars predicted to be observed by GLAST peaks at a lower 
7-ray luminosity towards the more populated region of radio pulsars than the distribution predicted 
for EGRET. Due to the greater sensitivity of GLAST, the group of radio quiet, 7-ray pulsars also 
moves toward lower 7-ray luminosities and more populated regions in the P — P space, but the 
radio-quiet pulsars are younger. 

The radio-quiet pulsar, Geminga, observed by EGRET is located in the GLAST region of 
radio-quiet pulsars in Fig. 10, but somewhat removed from the EGRET region. The results of 
the simulation seem to agree fairly well with the observations made by EGRET though, perhaps, 
predicting a few too many radio-loud, 7-ray pulsars. As mentioned before, many of the EGRET 
unidentified sources are expected to be radio-quiet, 7-ray pulsars. 

The results predicted for GLAST are interesting in that about as many radio-quiet, 7-ray 
pulsars are expected to be observed as radio-loud, 7-ray pulsars. Given the flux thresholds we have 
used, these are objects observed as point sources. In order to understand why GLAST is predicted 
to detect a larger ratio of radio-quiet to radio-loud, 7-ray pulsars than EGRET, we present in Figure 
11 the distribution of the distance of the pulsar from Earth for each of these groups. We simulated 
a group of 10,000 radio pulsars and normalized the distributions to the number (445) of observed 
radio pulsars in the selected group. As a result, one can see fractional numbers of pulsars in the 
distributions. The distribution of observed radio pulsars is presented as a shaded histogram. The 
predicted distributions of radio-quiet (thin) and radio-loud (thick), 7-ray pulsars for GLAST are 
displayed by the plain histograms. The predicted distributions of radio-quiet (diagonal pattern) and 
radio-loud (brick pattern), 7-ray pulsars are indicated for EGRET. The distribution of radio-loud, 
7-ray pulsars detected by both EGRET and GLAST have distributions with similar shapes as the 
observed radio distribution. Due to significantly different 7-ray flux thresholds, the distributions of 
radio-quiet, 7-ray pulsars are quite different for EGRET and GLAST. Thus the simulation suggests 
that due to the increased sensitivity, GLAST will be able to detect more pulsars that are further 
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away than even those detected by the radio surveys used in this study. This may change with 
the observations from the Parkes multibeam pulsar survey (Manchester et al. 2001) that has been 
finding more young, distant pulsars. 

In Table 8, we summarize the number of radio-quiet and radio-loud, 7-ray pulsars simulated 
for the two main cases explored in this study in which we have assumed no magnetic field decay 
requiring a pulsar death valley and field decay without a death valley. The results for the case of 
no field decay without a valley are also included for comparison. In addition, we have indicated the 
number of pulsed sources that GLAST would observe from the radio-quiet, 7-ray pulsar group. 

Table 8 



Simulated Pulsar Statistics 





EGRET 


GLAST 




Neutron Star Birth 


Case 


Radio Quiet 


Radio Loud 


Radio Quiet (Pulsed) 


Radio Loud 


Rate (per century) 


No decay 












with valley 


1 


7 


74 (7) 


76 


1 


No decay 












no valley 


1 


3 


35 (3) 


38 


0.5 


Decay 












no valley 


2 


9 


101 (9) 


90 


2 



These simulations are able to suggest a neutron star birth rate. As mentioned previously, we have 
assumed a 1 steradian solid angle of the radio beam, so there is a factor 47r more pulsars whose 
beams do not point in the direction of the Earth than those used in the simulation. Wc indicate 
in the tabic the estimated neutron star birth rates per century for each of the indicated cases. We 
have corrected the birth rates assuming a Gaussian beam to correct the overall detection by the 
ratio of the detection volume to that of the actual volume that increases the birth rate by a factor 
of 1.4 (Arzoumanian, Chernoff k. Cordes 2001). 

4. Discussion 

We have modeled the Galactic population of radio pulsars using a Monte Carlo simulation 
in order to calculate the expected numbers of 7-ray pulsars detectable by EGRET and GLAST. 
This paper has focused on the 7-ray luminosity predicted by the polar cap model of Daugherty 
& Harding (1996) and Zhang &: Harding (2000). Our simulation predicts about the same number 
of radio-loud 7-ray pulsars (7) compared to the number (8) detected by EGRET. We predict that 
GLAST should detect on the order of bf 75-100 radio-loud pulsars, however this number could also 
change with a detailed treatment of geometry. A very interesting, and somewhat unexpected, result 
of our simulation is the prediction that GLAST will detect about the same number (74) of radio- 
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quiet as radio- loud pulsars (76), at least as point sources. Since we have assumed that the radio 
and 7-ray beams are identical in this study, this result implies that GLAST will be more sensitive 
than radio surveys with regard to the detection of young pulsars in the Galaxy. Even with the 
considerable model uncertainties, we can conclude that polar cap models predict a much smaller 
ratio of radio-quiet to radio-loud 7-ray pulsars than do outer gap model studies. For example, 
Zhang, Zhang &; Cheng (2000) predict that GLAST will detect only about 80 radio-loud 7-ray 
pulsars but about 1100 radio-quiet 7-ray pulsars. This ratio then, regardless of the exact numbers, 
will be an important discriminator between polar cap and outer gap models. 

We have found that magnetic field decay could have a significant effect on pulsar evolution 
and the formation of the observed pulsar P — P distribution. As mentioned previously, there are 
distinct differences between the observed distribution of Figure 4a and the simulation without field 
decay of Figure 4b. There is an over abundance of simulated pulsars with high fields and high 
periods. The apparent drop in observed pulsars in this region could be conceived as evidence for 
the decay of the magnetic field. Our preliminary calculations suggest that a decay constant of the 
order of 10^ years is required to effect such a distribution, as also suggested by the ages of these 
pulsars in this region. A recent study of the radio pulsar distribution by Tauris & Konar (2001) 
also concludes that some torque decay is required. Most studies of field decay estimate longer decay 
constants in fields below B ~ 10^^ G. However, we believe that this issue warrants further study. 

We can analytically calculate the time it takes a pulsar from birth with a given magnetic field 
to reach the multipole death line. Since we are assuming a constant birth rate, dividing this time by 
the maximum age of 10^ years used in the simulation provides a survival fraction of initial pulsars 
found to the left of the multipole line and is given by the expression 

Fs = 9.0 X 10-^B~2^'^. (42) 

From this equation, one can see that as the magnetic field increases by an order of magnitude, the 
survival fraction will decrease by an order of magnitude. This explains why the factor ^3 in Table 1, 
though much smaller than Ai , makes a significant contribution. The initial population is weighted 
by this primary distribution of the magnetic field. Looking at the wings of the field distribution say 
at 10"^^ G and 10^^ G, assuming no selection effects, these pulsars would be randomly distributed 
over time with 90% of them lying between the age lines of 10^ and 10^ years. However, the 
radio luminosity, being proportional to P^/^P~^, is dropping toward the multipole line. Therefore, 
pulsars with shorter periods to the left of the multipole line are more radio luminous and are 
seen scattered to the left of the line for an age of 10^ years. The simulation indicates a broader 
distribution in period along the B = 10^^ G line than the distribution along the B = 10^^ G line. 
The observed distribution shows the opposite effect in these regions with a narrower distribution 
along the B = 10^^ G line with very few pulsars with periods less than 0.2 seconds and a broader 
distribution along the line B = 10^^ G that cuts off the number of pulsars dramatically beyond a 2 
second period. The paucity of pulsars at low periods is maybe related to several effects that we have 
not accounted for in this simulation. Perhaps we have underestimated the actual minimum radio 
flux threshold given theoretically by equation (27) and plotted in Figure 1. The flux threshold 
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tends to increase dramatically for some surveys for periods less than 0.1 s (as in the cases the 
Molonglo 2 and Green Bank 2 surveys with rather large sampling rates) , which is where the observed 
distribution in Figure 4a begins to indicate fewer pulsars than the simulated distribution in Figure 
4b. This may be a result of the fact that we have left out the fifth term in equation (29), which 
becomes significant for low period pulsars and is more significant for the Green Bank 3 and Arecibo 
2 surveys. We hope to add this term to our next simulations when we also take into account the 
emission geometry. However, this is not a problem for the newer surveys with smaller sampling 
rates. With the exception of the Jodrell Bank 2 and Parkes 1 survey, the surveys were performed 
at low frequencies around 420 MHz. Young pulsars with short periods and large period derivatives 
are distant and in the Galactic plane where the effects of scattering and background temperatures 
require a more careful treatment than we have done in this simulation. 

On the other hand, the group of high field, large period pulsars missing in the observed 
distribution may be a manifestation of the effects of the geometry of the radio beam. Various 
studies (Arzoumanian, Chernoff &; Cordes 2001, Kijak & Gil 1998, Gil k. Han 1996 and Rankin 
1993) suggest that the pulse width of both the core and conal emission beams is proportional to 
p-i/2 Xhcrcforc, as the period increases, the solid angle of the emission decreases. The observed 
distribution in Figure 4a seems to indicate a decrease of pulsars above a period of about 2 seconds 
across magnetic field strengths higher than 10^^ G to the left of the multipole death line, which is 
also reflected in Figure 7 where the simulation predicts too many older pulsars. We expect that 
our forthcoming model with more realistic geometries will yield better agreement in this region. 

In our model, we introduced a death valley to drastically reduce the number of pulsars between 
the dipole death lines and the multipole death line. As the pulsars age along constant field lines 
between these lines both the 7-ray and radio luminosities are decreasing as well as the spin-down 
energy. Cascade simulations (Baring & Harding 2001) suggest that the density of electron-positron 
pairs is also decreasing in a similar fashion. While there has been no firm theoretical development 
to connect the pair density and the radio luminosity, one can speculate that there must be a definite 
relationship, and the radio luminosity must decrease in some manner with decreasing pair density. 
Perhaps the need for the death valley also reflects the geometric relationship of the radio beam to 
the period, as discussed above. 

There is also an interdependence of the derived velocity distribution, radio luminosity function 
and field evolution. For example, the velocity distribution predicted by Lyne & Lorimer (1994) does 
not fit the observed z distribution assuming the luminosity model of Narayan & Ostriker (1990). 
Using their distribution, we find a width of the z distribution that is twice that of the observed 
distribution as shown in Figure 9. Our simulations suggest that most of the pulsars can be accounted 
for with a smaller mean velocity of 170 km/s, which we found to be essentially independent 
on whether the field decays or not. Clearly there are pulsars with space velocities of at 
least one thousand km/s. However, we did not attempt to introduce a second distribution of high 
velocities as in the case of the study by Cordes & Chernoff (1998) and Arzoumanian, Chernoff & 
Cordes (2001). We have also been careful to check that the distance distribution of the simulated 
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pulsars agrees with the observed distribution as noted in Figure 7. If the pulsars are too radio 
bright, more pulsars further away will be detected in the simulation causing the distribution to 
shift further away. Assuming the distance model is correct, we adjusted the Narayan &: Ostriker 
(1990) luminosity law based on the observed luminosity distributions shown in Figure 8. The 
luminosity model of Narayan k, Ostriker (1990) is also somewhat dependent on their assumption of 
field decay, for which they studied several cases with decay constants ranging from 8 million to 11 
million years. Bhattacharya et al. (1992), using the same luminosity model, found similar results 
for decay constants of 10 and 100 million years, although the study of Hartman et al. (1997), which 
did not include magnetic field decay, also used this same luminosity law. So the luminosity model, 
velocity distribution and field evolution remain uncertain, especially without the inclusion of the 
geometry of the radio beams. 

A population of radio-quiet, 7-ray pulsars results very naturally from the different dependence 
upon the pulsar period and period derivative of the radio and 7-ray luminosities. With the same 
beam geometry for both radio and 7-ray emission, the model predicts a few radio-quiet, 7-ray 
pulsars, like Geminga, that EGRET should have observed as well as about the right number of 
radio-loud, 7-ray pulsars. The number if radio-quiet 7-ray pulsars is likely to increase with a 
more realistic treatment of beaming geometry. Within the EGRET Third Catalog, there are 171 
unidentified 7-ray point sources, many of which could be radio-quiet, 7-ray pulsars, an issue that 
will be settled by GLAST with its ability to perform period searches. 

If there were no selection effects imposed, the main concentration of pulsars would be near the 
multipole death line. Since the average radio luminosity is proportional to P^/^P~^, the observed 
radio population is pushed away from the multipole death line towards the upper left portion of 
the P — P diagram. The average 7-ray luminosity is proportional to P^/^P^V^ ^nd, therefore, 
drops faster than the radio luminosity pushing the population of 7-ray pulsars even higher toward 
the upper left, especially in the case of EGRET with higher thresholds than GLAST as seen in 
Figure 10. With greater sensitivities of GLAST, the 7-ray pulsar population moves more to the 
lower right to correspond more with the region of the population of radio pulsars as also noted in 
Figure lOd. However, with the lower 7-ray flux thresholds expected for GLAST, it will be able to 
observe more pulsars at larger distances than those observed by the radio telescopes that were used 
in the observations of the select group of pulsars in the catalog. As a result, the expected number 
of radio-quiet, 7-ray pulsars increases relative to the radio-loud, 7-ray pulsars from those observed 
by EGRET. 

The ratio of radio-quiet to radio-loud, 7-ray pulsars is strongly sensitive to the assumptions 
for the geometries of the radio and 7-ray beams. In this study, we assumed a simple geometry 
that accords with the polar cap model where the region of radio and 7-ray emission are tied to the 
magnetic polar cap area and are very similar. The outer gap model, on the other hand, suggests a 
very different geometry for the 7-ray beam, which originates on field lines of the opposite pole from 
that of the visible radio beam. These outer gap models (Romani & Yadigaroglu 1995, Cheng & 
Zhang 1998, Zhang, Zhang k, Cheng 2000) predict that GLAST will detect many more radio-quiet, 
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7-ray pulsars. In a model independent study, MacLaughlin & Cordes (2001) simply assume a large 
solid angle (27r) for the 7-ray beam and predict that GLAST will detect around 120 radio-loud 
and 750 radio-quiet pulsars. Polar cap models predict that older pulsars with ages greater than 
10^ years and with longer periods should emit 7 rays, while the outer gap models predict that 
such pulsars will be radio-quiet (Chen k, Ruderman 1993). bf While the simple model studied here 
provides strong diagnostics, in the near future, we hope to perform simulations with a model that 
includes a more realistic geometry for the radio as well as the 7-ray beams and should provide 
further understanding that may help differentiate between pulsar models. 
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Fig. 1. Radio flux thresholds, Smim at the observing frequency for each of the eight surveys used 
in the simulation as a function of the pulsar period assuming a sky temperature of 150 K at 408 
MHz and a dispersion measure of 200 cm~^ • pc. 
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Fig. 2. — The average 7-ray energy (contours in GeV) from equation (39) is plotted (a) as a function 
of spectral index and the high energy cutoff assuming a threshold energy, eth = 100 MeV, typical 
for EGRET, bf (b) The logarithm of the average energy (contours) in GcV as a function of the 
logarithms of the period derivative and the period from equations (39) to (41). 




Fig. 3. — Aitoff plots of the observed pulsars (a) and of the simulated pulsars (b) for the case of 
no field decay and a death valley. 
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Period Period 

Fig. 4. — Distributions of observed pulsars (a) and simulated pulsars (b) for the case of no field 
decay and a death valley as a function of the period derivative and period (in seconds) of the 
pulsars. Solid dots indicate radio pulsars. Solid triangles represent radio-loud, 7-ray pulsars and 
open circles symbolize radio-quiet, 7-ray pulsars observed (a) and predicted (b) for EGRET. 
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Fig. 5. — Distributions of simulated pulsars as a function of the period derivative and period 
for case of no field decay and no death valley. Solid dots indicate radio pulsars. Solid triangles 
represent radio-loud, 7-ray pulsars and open circles symbolize radio-quiet, 7-ray pulsars predicted 
for EGRET. 
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Fig. 6. — Distributions of simulated pulsars as a function of the period derivative and period of 
simulated pulsars for case of magnetic field decay, assuming a time constant of 5 x 10^ years, and 
no pulsar death valley. Solid dots indicate radio pulsars. Solid triangles represent radio-loud, 7-ray 
pulsars and open circles symbolize radio-quiet, 7-ray pulsars predicted for EGRET. 
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Fig. 7. — Distributions of various pulsar characteristics indicated as shaded histograms (observed 
pulsars) and plain histograms (simulated pulsars). Thick histograms represent the distributions of 
the simulated pulsars assuming no field decay and pulsar death valley, while the thin histograms 
result from the case assuming field decay and no death valley. 
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Fig. 8. — The logarithm of the radio luminosity as a function of the log(i^/P^) for the cases of an 
intercept of 7.2 (a) and 6.64 (c). Pulsars are simulated assuming no magnetic field decay and the 
presence of a death valley between the death lines. Observed and simulated pulsars represented by 
solid dots and open triangles, respectively. Linear fits of the observed and simulated distributions 
are represented by solid and dotted lines, respectively The distributions of the logarithm of the 
radio luminosity at 400 MHz for observed (shaded histogram) and predicted (plain histograms) 
pulsars for the cases of an intercept of 7.2 (b) and 6.64 (d). We find that similar results are 
obtained assuming field decay. 



-36- 



Q. 



CD 

E 



80 R 



60 - 



40 - 



20 - 



Observed 

Calculated with C=V/120 
- Calculated with i;=V/350 




Fig. 9. — The z distribution from the Galactic disk for observed (shaded histograms) and simulated 
pulsars for the cases of ^ = ■u/120 (thick line) and of ^ = v/350 (thin line) assuming no field decay 
and a death valley. For the case of field decay, the simulated distribution is slighty 
narrower with a width of 0.14 it 0.01 kpc. 
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Fig. 10. — Distribution of radio-quiet (upper panels a and b) and radio- loud (lower panels c and d), 
7-ray pulsars as a function of period derivative and period for EGRET and GLAST. Observed and 
simulated pulsars are represented by solid triangles and solid circles, respectively. The simulated 
distributions are for the case of no magnetic field decay and the presence of a pulsar death valley. 
Contours represent the distributions for a simulation of 10,000 radio pulsars in order to obtain 
smoother distributions. The background gray scale plots represent the average 7-ray (upper panels 
a and b) and radio (lower panels c and d) luminosities. 
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Fig. 11. — Distributions of the distance from Earth of observed radio pulsars (shaded histograms), 
radio-loud (thick histograms) and radio-quiet (thin histograms), gamma pulsars for GLAST), radio- 
loud (brick patterned histograms) and radio-quiet (diagonally patterned histograms), gamma pul- 
sars for EGRET. The simulated distributions are for the case of no magnetic field decay and the 
presence of a pulsar death valley. 



